In [1]:
%matplotlib inline
In [136]:
#############
# Count number of strings in list that are the same
#
#
#############
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]
from collections import defaultdict
import os
import numpy as np
import re
#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')
#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]
for i, c in enumerate(station_data):
if c[3]:
match_station= re.compile(r"%s" % c[3][0][0][0:5])
if match_station.search(str(station_list_pick)) is not None:
#print st_lat_c
st_lat_c.append(c[1])
st_lon_c.append(c[2])
st_nam_c.append(c[0].lower().title())
#st_cnt_c[i]=str()
#st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
#st_namcnt_c.append(st_nam_c[i])
#print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances
slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset
lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]
initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\
m = Basemap(projection='cyl',\
llcrnrlat=-10.,urcrnrlat=40.,\
llcrnrlon=60.,urcrnrlon=100.,\
rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
for i,j,s in zip(x, y, st_nam_c):
plt.text(i, j, '$%s$' % s, fontsize=14)
plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png', format='png', bbox_inches='tight')
m.plot(x, y_line)
plt.show()
for x in sorted(st_nam_c):
print x
In [137]:
#############
# Count number of strings in list that are the same
#
#
#############
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]
station_list_pick=[43371, 43353, 43285, 43192, 43003]
# Thiruvanthanapuram, Cochin, Panambur, Panjim, Bombay
from collections import defaultdict
import os
import numpy as np
import re
#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')
#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]
for i, c in enumerate(station_data):
if c[3]:
match_station= re.compile(r"%s" % c[3][0][0][0:5])
if match_station.search(str(station_list_pick)) is not None:
#print st_lat_c
st_lat_c.append(c[1])
st_lon_c.append(c[2])
st_nam_c.append(c[0].lower().title())
#st_cnt_c[i]=str()
#st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
#st_namcnt_c.append(st_nam_c[i])
#print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances
slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset
lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]
initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\
m = Basemap(projection='cyl',\
llcrnrlat=-10.,urcrnrlat=40.,\
llcrnrlon=60.,urcrnrlon=100.,\
rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
for i,j,s in zip(x, y, st_nam_c):
plt.text(i, j, '$%s$' % s, fontsize=14)
plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png', format='png', bbox_inches='tight')
m.plot(x, y_line)
plt.show()
for x in sorted(st_nam_c):
print x
In [138]:
#############
# Count number of strings in list that are the same
#
#
#############
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]
station_list_pick=[43003, 43014, 42867]
# Bombay, Aurangabd, Nagpur
from collections import defaultdict
import os
import numpy as np
import re
#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')
#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]
for i, c in enumerate(station_data):
if c[3]:
match_station= re.compile(r"%s" % c[3][0][0][0:5])
if match_station.search(str(station_list_pick)) is not None:
#print st_lat_c
st_lat_c.append(c[1])
st_lon_c.append(c[2])
st_nam_c.append(c[0].lower().title())
#st_cnt_c[i]=str()
#st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
#st_namcnt_c.append(st_nam_c[i])
#print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances
slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset
lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]
initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\
m = Basemap(projection='cyl',\
llcrnrlat=-10.,urcrnrlat=40.,\
llcrnrlon=60.,urcrnrlon=100.,\
rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
for i,j,s in zip(x, y, st_nam_c):
plt.text(i, j, '$%s$' % s, fontsize=14)
plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png', format='png', bbox_inches='tight')
m.plot(x, y_line)
plt.show()
for x in sorted(st_nam_c):
print x
In [122]:
def calculate_initial_compass_bearing(pointA, pointB):
"""
Calculates the bearing between two points.
The formulae used is the following:
θ = atan2(sin(Δlong).cos(lat2),
cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong))
:Parameters:
- `pointA: The tuple representing the latitude/longitude for the
first point. Latitude and longitude must be in decimal degrees
- `pointB: The tuple representing the latitude/longitude for the
second point. Latitude and longitude must be in decimal degrees
:Returns:
The bearing in degrees
:Returns Type:
float
"""
if (type(pointA) != tuple) or (type(pointB) != tuple):
raise TypeError("Only tuples are supported as arguments")
lat1 = math.radians(pointA[0])
lat2 = math.radians(pointB[0])
diffLong = math.radians(pointB[1] - pointA[1])
y = math.sin(diffLong) * math.cos(lat2)
x = math.cos(lat1) * math.sin(lat2) - (math.sin(lat1)
* math.cos(lat2) * math.cos(diffLong))
initial_bearing = math.atan2(y, x)
# Now we have the initial bearing but math.atan2 return values
# from -180° to + 180° which is not what we want for a compass bearing
# The solution is to normalize the initial bearing as shown below
initial_bearing = math.degrees(initial_bearing)
compass_bearing = (initial_bearing + 360) % 360
return initial_bearing, compass_bearing
def compass_to_maths_bearing(initial_bearing):
'''
Converts degrees positive clockwise from north (in range +/- 180, or 0-360) bearing to
positive anticlockwise from east maths bearing
'''
if initial_bearing > 180:
initial_bearing=initial_bearing-360
maths_bearing=(-initial_bearing)+90
return maths_bearing
In [132]:
import numpy as np
from math import *
import math
slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset
lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]
initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
#u_rot = U * cos(brng) + V * sin(brng)
#v_rot = -U * sin(brng) + V * cos(brng)
Out[132]:
θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )
Bearing
In general, your current heading will vary as you follow a great circle path (orthodrome); the final heading will differ from the initial heading by varying degrees according to distance and latitude (if you were to go from say 35°N,45°E (≈ Baghdad) to 35°N,135°E (≈ Osaka), you would start on a heading of 60° and end up on a heading of 120°!).
This formula is for the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a great-circle arc will take you from the start point to the end point:1 Formula: θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ ) JavaScript: (all angles in radians)
var y = Math.sin(λ2-λ1) Math.cos(φ2); var x = Math.cos(φ1)Math.sin(φ2) - Math.sin(φ1)Math.cos(φ2)Math.cos(λ2-λ1); var brng = Math.atan2(y, x).toDegrees();
Excel: (all angles in radians) =ATAN2(COS(lat1)SIN(lat2)-SIN(lat1)COS(lat2)COS(lon2-lon1), SIN(lon2-lon1)COS(lat2)) *note that Excel reverses the arguments to ATAN2 – see notes below
Since atan2 returns values in the range -π ... +π (that is, -180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with −ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is (floating point) modulo.
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
N.B - lats and lons have to be converted to radians first